class: center, middle, inverse, title-slide # Taller de R: Estadística y programación ## Lectura 13: Datos raster ### Eduard Martínez ### Universidad de los Andes |
ECON-1302
--- name: contenido # Ya vimos.. - [Operaciones geométricas con objetos sf](https://github.com/taller-R/clase-12) -- # Hoy veremos 1. [Introducción a datos raster](#intro) 2. [Importar datos raster](#read) 3. [Operaciones geométricas con raster](#ope) 4. [Raster de varios layers](#layers) <!-------------------------------> <!--- Introducción a datos raster ---> <!-------------------------------> --- class: inverse, center, middle name: intro # [1. Introducción a datos raster](#contenido) <html><div style='float:left'></div><hr color='#FF007F' size=1px width=796px></html> --- # 1.1. Qué es un raster? - Un raster es una matriz de celdas (o píxeles) organizadas en filas y columnas (o una cuadrícula) en la que cada celda contiene un valor que representa información. <div align="center"> <img src="pics/raster.png" height=400> </div> Tomado de: [https://www.neonscience.org](https://www.neonscience.org/resources/learning-hub/tutorials/dc-raster-data-r) --- # 1.2. Elementos de un raster * Dimensión o numero de filas y columnas, así como el numero de bandas que tiene el raster * La resolución de cada pixel * Tipo de archivo * CRS * Información dentro de los pixelex # 1.3. Extenciones de un raster * .tif * .nc * .nc4 --- # 1.4. Algunas fuentes de información de datos raster - **Luces nocturnas mensuales 2012-2019** Pueden descargarse en [https://ngdc.noaa.gov](https://ngdc.noaa.gov/eog/viirs/download_dnb_composites.html) o en [https://payneinstitute.mines.edu](https://payneinstitute.mines.edu/eog/) - **Deforestación (SIAC)** Pueden descargarse en [http://www.siac.gov.co](http://www.siac.gov.co) - **Deforestación (Hansen)** Pueden descargarse en [https://data.globalforestwatch.org](https://data.globalforestwatch.org/datasets/14228e6347c44f5691572169e9e107ad) - **GES-DISC** Una de las principales fuentes de datos raster a nivel mundial es *EARTH DATA* de la NASA y puedes descargar los raster en [https://disc.gsfc.nasa.gov](https://disc.gsfc.nasa.gov). - **Librería getSpatialData** La librería getSpatialData provee raster de los proyetos Sentinel-1, Sentinel-2, Sentinel-3, Sentinel-5P, Landsat 8 OLI, Landsat ETM, Landsat TM, Landsat MSS, MODIS (Terra & Aqua) y SRTM DEMs de la NASA. Puede consultar el repositorio de la libreria en [https://github.com/16EAGLE/getSpatialData](https://github.com/16EAGLE/getSpatialData). <!-------------------------------> <!--- Importar datos raster ---> <!-------------------------------> --- class: inverse, center, middle name: read # [2. Importar datos raster](#contenido) <html><div style='float:left'></div><hr color='#FF007F' size=1px width=796px></html> --- # 2.1. Cargemos los datos (...) ```r # llamar las librerias require(raster) # Leer raster deforestacion = raster("data/input/siac/magdalena_deforestacion_1990_2000.tif") ``` ``` ## Warning in showSRID(uprojargs, format = "PROJ", multiline = "NO", prefer_proj ## = prefer_proj): Discarded datum Unknown based on GRS80 ellipsoid in Proj4 ## definition ``` ```r # atributos del raster deforestacion ``` ``` ## class : RasterLayer ## dimensions : 2505, 2704, 6773520 (nrow, ncol, ncell) ## resolution : 0.0002753217, 0.0002713948 (x, y) ## extent : -74.28793, -73.54346, 10.42546, 11.1053 (xmin, xmax, ymin, ymax) ## crs : +proj=longlat +ellps=GRS80 +no_defs ## source : /Users/eduard-martinez/Dropbox/teaching/Taller de R/GitHub/clases/clase-13/data/input/SIAC/magdalena_deforestacion_1990_2000.tif ## names : magdalena_deforestacion_1990_2000 ## values : 1, 5 (min, max) ``` --- # 2.2. Veamos el raster ```r # plot basico del raster plot(deforestacion) ``` <div align="center"> <img src="pics/Rplot.png" height=450> </div> --- # 2.3. Geometría del raster ```r # nombres de las bandas names(deforestacion) ``` ``` ## [1] "magdalena_deforestacion_1990_2000" ``` ```r names(deforestacion) = "cobertura_vegetal" names(deforestacion) ``` ``` ## [1] "cobertura_vegetal" ``` ```r # Extension st_bbox(deforestacion) ``` ``` ## xmin ymin xmax ymax ## -74.28793 10.42546 -73.54346 11.10530 ``` ```r # Proyeccion deforestacion@crs ``` ``` ## CRS arguments: +proj=longlat +ellps=GRS80 +no_defs ``` ```r crs(deforestacion) ``` ``` ## CRS arguments: +proj=longlat +ellps=GRS80 +no_defs ``` --- # 2.4. Valores del raster (...) ```r # Valores de los pixeles minValue(deforestacion$cobertura_vegetal) ``` ``` ## [1] 1 ``` ```r maxValue(deforestacion$cobertura_vegetal) ``` ``` ## [1] 5 ``` ```r values(deforestacion) %>% head() ``` ``` ## [1] NA NA NA NA NA NA ``` ```r # Cambiar los values de los pixeles values(deforestacion)[1] <- 0 values(deforestacion)[is.na(values(deforestacion))==T] <- 0 values(deforestacion) %>% head() ``` ``` ## [1] 0 0 0 0 0 0 ``` --- # 2.4. Valores del raster (...) ```r # Descriptivas values(deforestacion) %>% table() ``` ``` ## . ## 0 1 2 3 4 5 ## 76469 2432531 200575 88669 26717 3948559 ``` ```r values(deforestacion) %>% summary() ``` ``` ## Min. 1st Qu. Median Mean 3rd Qu. Max. ## 0.000 1.000 5.000 3.388 5.000 5.000 ``` ```r value_raster = values(deforestacion) str(value_raster) ``` ``` ## num [1:6773520] 0 0 0 0 0 0 0 0 0 0 ... ``` <!-------------------------------> <!--- Operaciones geométricas con raster ---> <!-------------------------------> --- class: inverse, center, middle name: ope # [3. Operaciones geométricas con raster](#contenido) <html><div style='float:left'></div><hr color='#FF007F' size=1px width=796px></html> --- # 3.1. Importar raster de luces ```r # leer raster luces = raster(x = 'data/input/NOAA/colombia_202004.tif') luces ``` ``` ## class : RasterLayer ## dimensions : 2989, 2918, 8721902 (nrow, ncol, ncell) ## resolution : 0.004166667, 0.004166667 (x, y) ## extent : -79.00625, -66.84792, 0.002082733, 12.45625 (xmin, xmax, ymin, ymax) ## crs : +proj=longlat +datum=WGS84 +no_defs ## source : /Users/eduard-martinez/Dropbox/teaching/Taller de R/GitHub/clases/clase-13/data/input/NOAA/colombia_202004.tif ## names : colombia_202004 ## values : -0.43, 1389.84 (min, max) ``` ```r # leer polygono cataca = st_read(dsn= 'data/input/mgn/MGN_Municipio.shp') %>% subset(MPIO_CCDGO %in% c('053')) ``` ``` ## Reading layer `MGN_Municipio' from data source `/Users/eduard-martinez/Dropbox/teaching/Taller de R/GitHub/clases/clase-13/data/input/mgn/MGN_Municipio.shp' using driver `ESRI Shapefile' ## Simple feature collection with 30 features and 12 fields ## geometry type: POLYGON ## dimension: XY ## bbox: xmin: -74.9471 ymin: 8.912502 xmax: -73.54336 ymax: 11.34912 ## geographic CRS: WGS 84 ``` ```r cataca ``` ``` ## Simple feature collection with 1 feature and 12 fields ## geometry type: POLYGON ## dimension: XY ## bbox: xmin: -74.25773 ymin: 10.42533 xmax: -73.54336 ymax: 10.8791 ## geographic CRS: WGS 84 ## OBJECTID DPTO_DPTO_ MPIO_CCDGO MPIO_CNMBR ## 3 659 47 053 ARACATACA ## MPIO_CRSLC MPIO_NANO_ MPIO_NAREA MPIO_CSMBL ## 3 Ordenanza 47 del 28 de Abril de 1915 1993 1746.793 4 ## MPIO_CCNCT MPIO_NANO SHAPE_Leng SHAPE_Area geometry ## 3 47053 2012 2.650507 0.1443276 POLYGON ((-73.72812 10.8626... ``` --- # 3.2. Cliping un raster(...) ```r l_cataca = crop(luces,cataca) # Qué debo revisar primero? ggplot() + geom_tile(data = as.data.frame(l_cataca, xy=TRUE), aes(y=y,x=x,fill=colombia_202004)) + scale_fill_distiller(palette='Spectral',na.value = 'gray') + geom_sf(data =cataca,color = 'black',fill=NA) + theme_bw() ``` <img src="data:image/png;base64,#intro-to-class_files/figure-html/unnamed-chunk-7-1.png" style="display: block; margin: auto;" /> --- # 3.2. Cliping un raster(...) ```r l_cataca = crop(luces,cataca) %>% mask(cataca) ggplot() + geom_tile(data = as.data.frame(l_cataca, xy=TRUE), aes(y=y,x=x,fill=colombia_202004)) + scale_fill_distiller(palette='Spectral',na.value = 'gray') + geom_sf(data =cataca,color = 'black',fill=NA) + theme_bw() ``` <img src="data:image/png;base64,#intro-to-class_files/figure-html/unnamed-chunk-8-1.png" style="display: block; margin: auto;" /> --- # 3.3. Extraer los valores del raster ```r values(l_cataca) %>% .[is.na(.)==F] %>% head() ``` ``` ## [1] 0.00 0.00 0.00 0.00 0.21 0.00 ``` ```r data = values(l_cataca) %>% .[is.na(.)==F] summary(data) ``` ``` ## Min. 1st Qu. Median Mean 3rd Qu. Max. ## -0.3600 0.2600 0.3300 0.3862 0.4000 19.7600 ``` --- # 3.4. Raster a datos vectoriales (...) ```r # Raster a puntos point = rasterToPoints(l_cataca,spatial = T) %>% st_as_sf() class(point) ``` ``` ## [1] "sf" "data.frame" ``` ```r ggplot() + geom_sf(data = point , aes(color=colombia_202004),size=0.1) ``` <img src="data:image/png;base64,#intro-to-class_files/figure-html/unnamed-chunk-10-1.png" style="display: block; margin: auto;" /> --- # 3.4. Raster a datos vectoriales (...) ```r # Raster a polygonos polygon = rasterToPolygons(l_cataca) %>% st_as_sf() ggplot() + geom_sf(data = polygon , aes(fill=colombia_202004),size=0) + scale_fill_distiller(palette='Spectral',na.value = 'gray') ``` <img src="data:image/png;base64,#intro-to-class_files/figure-html/unnamed-chunk-11-1.png" style="display: block; margin: auto;" /> <!-------------------------------> <!--- Raster de varios layers ---> <!-------------------------------> --- class: inverse, center, middle name: layers # [4. Raster de varios layers](#contenido) <html><div style='float:left'></div><hr color='#FF007F' size=1px width=796px></html> --- # 4.1. Cómo se ve un raster RGB? <div align="center"> <img src="pics/rgb_raster.png" height=500> </div> Tomado de: [https://www.neonscience.org](https://www.neonscience.org) --- # 4.2. Importar raster ```r GDALinfo("data/input/neon/HARV_RGB_Ortho.tif") ``` ``` ## rows 2317 ## columns 3073 ## bands 3 ## lower left origin.x 731998.5 ## lower left origin.y 4712956 ## res.x 0.25 ## res.y 0.25 ## ysign -1 ## oblique.x 0 ## oblique.y 0 ## driver GTiff ## projection +proj=utm +zone=18 +datum=WGS84 +units=m +no_defs ## file data/input/neon/HARV_RGB_Ortho.tif ## apparent band summary: ## GDType hasNoDataValue NoDataValue blockSize1 blockSize2 ## 1 Float64 TRUE -1.7e+308 1 3073 ## 2 Float64 TRUE -1.7e+308 1 3073 ## 3 Float64 TRUE -1.7e+308 1 3073 ## apparent band statistics: ## Bmin Bmax Bmean Bsd ## 1 0 255 NaN NaN ## 2 0 255 NaN NaN ## 3 0 255 NaN NaN ## Metadata: ## AREA_OR_POINT=Area ``` --- # 4.3. Veamos que tenemos ```r banda_r <- raster(x = "data/input/neon/HARV_RGB_Ortho.tif") # R carga por 'default' la banda 1, es decir la roja paleta_col <- gray.colors(n = 100, start = 0.0,end = 1.0,alpha = NULL) plot(banda_r, col=paleta_col, axes=FALSE, main="Imagen RGB - Banda 1 (roja)") ``` <img src="data:image/png;base64,#intro-to-class_files/figure-html/unnamed-chunk-13-1.png" style="display: block; margin: auto;" /> --- # 4.4. Values de la banda En un RGB podemos tener 255*255*255 posibles combinaciones, es decir 16.581.375 colores. ```r minValue(banda_r) ``` ``` ## [1] 0 ``` ```r maxValue(banda_r) ``` ``` ## [1] 255 ``` ```r values(banda_r) %>% summary() ``` ``` ## Min. 1st Qu. Median Mean 3rd Qu. Max. ## 0.0 81.0 108.0 102.2 128.0 255.0 ``` --- # 4.5. cargar las 3 bandas ```r RGB_stack = stack("data/input/neon/HARV_RGB_Ortho.tif") plotRGB(RGB_stack, r = 1, g = 2, b = 3) ``` <img src="data:image/png;base64,#intro-to-class_files/figure-html/unnamed-chunk-15-1.png" style="display: block; margin: auto;" /> --- # 4.6. Inspeccionemos el objeto (...) ```r # tipo de objeto class(RGB_stack) ``` ``` ## [1] "RasterStack" ## attr(,"package") ## [1] "raster" ``` ```r dim(RGB_stack) ``` ``` ## [1] 2317 3073 3 ``` ```r # nombres de las bandas names(RGB_stack) ``` ``` ## [1] "HARV_RGB_Ortho.1" "HARV_RGB_Ortho.2" "HARV_RGB_Ortho.3" ``` ```r names(RGB_stack) = c('red','green','blue') names(RGB_stack) ``` ``` ## [1] "red" "green" "blue" ``` --- # 4.6. Inspeccionemos el objeto (...) ```r RGB_stack@layers # ver las bandas ``` ``` ## [[1]] ## class : RasterLayer ## band : 1 (of 3 bands) ## dimensions : 2317, 3073, 7120141 (nrow, ncol, ncell) ## resolution : 0.25, 0.25 (x, y) ## extent : 731998.5, 732766.8, 4712956, 4713536 (xmin, xmax, ymin, ymax) ## crs : +proj=utm +zone=18 +datum=WGS84 +units=m +no_defs ## source : /Users/eduard-martinez/Dropbox/teaching/Taller de R/GitHub/clases/clase-13/data/input/NEON/HARV_RGB_Ortho.tif ## names : red ## values : 0, 255 (min, max) ## ## ## [[2]] ## class : RasterLayer ## band : 2 (of 3 bands) ## dimensions : 2317, 3073, 7120141 (nrow, ncol, ncell) ## resolution : 0.25, 0.25 (x, y) ## extent : 731998.5, 732766.8, 4712956, 4713536 (xmin, xmax, ymin, ymax) ## crs : +proj=utm +zone=18 +datum=WGS84 +units=m +no_defs ## source : /Users/eduard-martinez/Dropbox/teaching/Taller de R/GitHub/clases/clase-13/data/input/NEON/HARV_RGB_Ortho.tif ## names : green ## values : 0, 255 (min, max) ## ## ## [[3]] ## class : RasterLayer ## band : 3 (of 3 bands) ## dimensions : 2317, 3073, 7120141 (nrow, ncol, ncell) ## resolution : 0.25, 0.25 (x, y) ## extent : 731998.5, 732766.8, 4712956, 4713536 (xmin, xmax, ymin, ymax) ## crs : +proj=utm +zone=18 +datum=WGS84 +units=m +no_defs ## source : /Users/eduard-martinez/Dropbox/teaching/Taller de R/GitHub/clases/clase-13/data/input/NEON/HARV_RGB_Ortho.tif ## names : blue ## values : 0, 255 (min, max) ``` --- ### 4.7. Extraer atributos ```r point_RGB = rasterToPoints(RGB_stack,spatial = T) %>% st_as_sf() point_RGB ``` ``` ## Simple feature collection with 7120141 features and 3 fields ## geometry type: POINT ## dimension: XY ## bbox: xmin: 731998.6 ymin: 4712956 xmax: 732766.6 ymax: 4713535 ## CRS: +proj=utm +zone=18 +datum=WGS84 +units=m +no_defs ## First 10 features: ## red green blue geometry ## 1 0 1 0 POINT (731998.6 4713535) ## 2 2 0 10 POINT (731998.9 4713535) ## 3 6 9 1 POINT (731999.1 4713535) ## 4 0 0 0 POINT (731999.4 4713535) ## 5 16 5 17 POINT (731999.6 4713535) ## 6 0 0 0 POINT (731999.9 4713535) ## 7 0 4 4 POINT (732000.1 4713535) ## 8 6 2 0 POINT (732000.4 4713535) ## 9 1 1 0 POINT (732000.6 4713535) ## 10 5 0 7 POINT (732000.9 4713535) ``` <!---------------------> <!--- Hoy vimos ---> <!---------------------> --- class: inverse, center, middle name: view # [Hoy vimos...](#contenido1) <html><div style='float:left'></div><hr color='#FF007F' size=1px width=796px></html> --- # Hoy vimos... - ☑ [Introducción a datos raster](#intro) - ☑ [Importar datos raster](#read) - ☑ [Operaciones geométricas con raster](#ope) - ☑ [Raster de varios layers](#layers)